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Abstract 

Thermoacoustic computed tomography (thermoacoustic CT) has 
. the potential to become a mayor non-invasive medical imaging method. 

' In this paper we derive a general mathematical framework of a novel 

. measuring setup introduced in [P. Burgholzer, C. Hofcr. G. Paltauf. 

M. Haltmeier, and O. Scherzer, Thermoacoustic tomography with inte- 
grating area and line detectors, IEEE Transactions on Ultrasonics, Fer- 
roelectrics, and Frequency Control, 52 (2005)], that uses line shaped 
\Q . detectors instead of the usual point like ones. We show that the three 

' dimensional thermoacoustic imaging problem reduces to the mathc- 

■ matical problem of reconstructing the initial data of the two dimen- 

sional wave equation from boundary measurements of its solution. We 

S derive and analyze an analytic reconstruction formula which allows for 
fast numerical implementation. 
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1 Introduction and Motivation 

Thermoacoustic CT (also called opto- or photoacoustic CT) is a novel hybrid 
non-invasive imaging method with applications in different areas, e.g. in 
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medical diagnostics [15] or imaging of small animals [24]. It is based on the 
excitation of acoustic waves inside an investigated object when exposed to 
non-ionizing electromagnetic radiation [10, 8] and combines the advantages 
of purely optical imaging (high contrast) with ultrasonic high resolution 
[25, 4]. " 

So far the generated thermoacoustic waves are measured with several ul- 
trasonic transducers located outside the illuminated sample. The measured 
output of the ultrasonic transducers has been identified with the restriction 
of the thermoacoustic pressure field to a surface enclosing the object. Based 
on this point-data approximation the absorption density function can be 
reconstructed by solving the problem of recovering a function from its mean 
values over spherical surfaces [7, 25] . When conventional piezoelectric ultra- 
sonic transducers [25, 15] are used to approximate point-data, the necessity 
of using small detectors with high bandwidth is technically hardly realizable 
[6, 20]. 

To obtain high resolution the size of the detectors has to be taken into 
account when modeling the corresponding forward operator. Exact recon- 
struction formulas incorporating the detector size have been derived for large 
planar detectors in spherical geometry [9] and line detectors with cylindrical 
circular recording geometry [5, 19]. 

In this paper we establish a mathematical foundation of thermoacous- 
tic CT using line detectors in general recording geometry. We show that 
the three dimensional absorption density function can be reconstructed by 
solving the mathematical problem of recovering the initial data of the two 
dimensional wave equation from its solution, measured by an array of line de- 
tectors, on a recording curve C Tec . Therefore, line detectors offer a reduced 
numerical complexity. Finally three dimensional image reconstruction is 
achieved by rotating the array around a single axis. 

In the case where the recording curve C rec is a line, the initial data of the 
two dimensional wave equation can be recovered by a Fourier reconstruction 
formula [13, 14, 19]. In this article we present a rigorous mathematical anal- 
ysis of this formula, taking into account that the restriction of the solution 
of the wave equation to a line is not necessarily absolutely integrable. 

The paper is organized as follows. In Section 2 we recall the basic for- 
mulas of thermoacoustic CT and introduce the concept of line detectors. 
Section 3 is devoted to the decomposition of the three dimensional forward 
operator into a system of two dimensional operators corresponding to the 
two dimensional wave equation. In Section 4 we deal the case of linear 
recording geometry and present the analysis of the Fourier reconstruction 
formula. 
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2 Mathematical modeling of thermoacoustic CT 
with line detectors 



The basic principle of thermoacoustic CT is the generation of acoustic waves 
within an object by illuminating it with non- ionizing pulsed electromagnetic 
energy. 

Assume that a short electromagnetic pulse is emitted into a weakly ab- 
sorbing medium at time t = 0. The absorbed energy induces thermal heat- 
ing, which causes thermoelastic expansion and thereby an initial pressure 
distribution. This actuates small vibrations in the medium, which results in 
sound waves. 

The induced pressure field is mathematically modeled by a function p : 
R 3 x [0,oo) — ► R, where p(x,t) represents the acoustic pressure at position 
x S R 3 and time t =: t/c > 0. Here c denotes the speed of sound which is 
assumed to be constant. In this case the time varying pressure field satisfies 
the three dimensional wave equation (see [8, 25, 10]) 

^- A x Jp(x,i) = 0, (x,i) GR 3 x (0,oo) (1) 
with initial conditions 

p(x,0)=/(x), xGR 3 , (2) 

^(x,0)=0, xGR 3 . (3) 

The mathematical task in thermoacoustic CT is to reconstruct the initial 
data / using data of the solution of the three dimensional wave equation 
(1), (2), (3) gathered with several detectors located outside the investigated 
sample. 

The measured output of a conventional ultrasonic transducer is the total 
pressure 



t4 J p(x,t)dS(x) 



acting on the surface S of the ultrasonic transducer convolved in time with 
the temporal impulse response function of the detector [11]. Here dS denotes 
the two dimensional surface measure and |<S| the surface area of S. In 
thermoacoustic CT the temporal impulse response function of the ultrasonic 
detector is assumed to be an approximation of the ^-distribution and the 
surface area \S\ to be small. In this case the measurement data can be 
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Illumination 




Figure 1: Conventional thermoacoustic scanning system. The non- 
uniform absorbtion causes outgoing acoustic waves that are measured with 
an ultrasonic transducer located at different positions. 

approximated by 

p(x rcc ,t) = i J p(x,t)dS(x) 

where x rec is the center of the surface S. 

In practice conventional piezoelectric ultrasonic transducers placed on 
a recording surface 5 rec are used to approximate point-like detectors [15, 
25]. In this case thermoacoustic CT leads to the mathematical problem 
of recovering a function from its integrals over spheres with, centers x rcc 
(locations of the point-like detectors) on S vec [25, 7, 14]. Every ultrasonic 
transducers has a finite size (typically in the range of mm) and therefore 
algorithms that are based on the assumption that point measurement data 
is collected give blurred reconstructions. 

Line detectors measure the total pressure 

J dL{x) 

on a line L [5, 19]. In practical experiments they can be realized by a 
thin laser beam that is either part of Fabry-Perot [5] or Mach-Zehnder [19] 
interferometer. Details on the technical realization of line detectors can be 
found in [5, 19]. We just note here that optical sensors offer high bandwidth 
and high signal to noise ratio. 
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Data acquisition 

Throughout this paper we assume that the initial data / is supported in the 
open set V rec . We call V rcc recording volume, if V rcc is invariant with respect 
to rotations around the e3 = (0, 0, 1) axis, i.e. 

Free = {-Rx : X G Free} 

holds for all rotations R G SO (3) with fix point set Me3. Important examples 
for recording volumes are the half space M 2 x (0, oo), the cylinder flxR, 
where D denotes the unit disc in R 2 . We set D vec := V rec H ({0} x R 2 ). 

The measured data is collected with an array of parallel line detectors 
that are tangential to V-eci orthogonal to e3 and rotated around the axis e3. 
The rotation around the axis e3 is parameterized by <x G S 1 and defined 
by 

i? CT ei = ei(er) := (01, a 2 , 0) , 
-Ro-e 2 = e 2 (cr) := (-cr 2 , 0i, 0) , 
Rcre 3 = e 3 . 

Here (ei,e2,e3) denotes the standard basis of M 3 and (ei(cr), e2(<x), e 3 ) is 
a positively oriented orthonormal basis. In practical applications the mea- 
surements, for fixed <x, are performed either with a array or by moving a 
single line detector along a recording curve C rec . 
Let 

L(a, y) := {sei(er) + yie 2 (er) + y 2 e 3 :s£l} 

be a line, where y = (2/1,2/2) G R 2 - The positions of the line detectors are 
given by L(<x,yre C ) with y rec G C TCC C dD rec which is called recording curve, 
see Figure 2. 

Since V-cc is rotationally invariant with respect to e3 and C rec C dD rcc 
all lines -L(<r,y rec ) (locations of the line detectors) are outside the support 
of /. The forward operator 

V : C °°(Vrec) -> C°° (S 1 X C r cc x (0, 00)) 

/ » V(f) := ( (<r, y rC c, t) ^ I p(x, t) dL(x) 

y J L((T,y rec ) 

collects the data measured by all line detectors L(a,y Tec ) for given initial 
data /. Here p(x, t) denotes the unique solution of (1), (2), (3). We note 
that (1), (2), (3) is a well posed problem and therefore V(f) can be cal- 
culated stable for arbitrary initial data /. Thermoacoustic CT with line 
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detectors deals with the solution of the operator equation V(f) = g with 
given (potentially noisy) measured data g. Therefore we study uniqueness, 
stability and explicit inversion of V . In this paper we focus on the derivation 
and implementation of inversion formulas. 



e 3 




► e2 



Figure 2: Scanning geometry for thermoacoustic CT using line de- 
tectors. The line detectors (orthogonal to the drawing plane) are uniformly 
distributed on the recording curve C rec and rotated around a single axis of 
rotation e3, indicated by parameter cr. 



3 Decomposition of the three dimensional forward 
operator into a system of two dimensional oper- 
ators 

Let Vr ec be a recording volume and C rec C dD rec a recording curve. In this 
Section we prove that the corresponding forward operator V decomposes into 
the product of the two dimensional Radon transform 1Z and the operator 

Q:C co (D rcc )^C oo (C rcc x(0,cx))) 

F^Q{F) := ((y rec ,t) ^ P(y«c,t)) 

that maps a planar function F € Cq^(D tcc ) onto the solution P of the two 
dimensional wave equation 

J^-A y )i> = 0, (y,t) £K 2 x (0,oo) (4) 
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with initial conditions 

P(y,0)=F(y), y G R 2 , (5) 
dP 

^(y,o)=o, yet 2 (6) 

restricted to the recording curve C rec . Here A y denotes the Laplacian in R 2 . 
We use the following notation: 

1. Define 

K® J 3 : C °°(y rec ) C^(S l x R 2 ) 

f»(K®Is)(f) := ((r,(r,z)) » (K(f z ))(r,rj) 

where 1Z denotes the classical two dimensional Radon transform [17] 
defined by 

K : C^(M. 2 ) -> C^iS 1 x R) 

<p i— > TZ((p) ■= mt, r) i— > y y?(sr + rr 1 ) ds^j 

and / 2 the function / for fixed third parameter 2, i.e. 

f z (x,y) := f(x,y,z) , with z e R fixed . 

Here r -1 " is a unit vector orthogonal to r. The notation indicates 
that the operator TZ®!^ operates on a function / defined on R 3 by 
applying the Radon transform to the first two components for fixed 
third component. 

2. For a fixed reS 1 and a function G £ Cq (5 1 x -D rcc ) define 

GV(xi,x 2 ) := G(r,x 1 ,x 2 ) ■ 

Furthermore, we define an operator I® Q by 

1®Q: C^iS 1 x £> rec ) C°° (S 1 x C rcc x (0, 00)) 

G^{1®Q){G) := ((r,y iec ,t)^ ( Q(G T ))(y rcc ,t)) 

The main feature of using line detectors is based on the following de- 
composition: 
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Theorem 1. Let f G C$°(V Iec ). Then 

V(f) = (I®Q)o(n®I 3 )(f). (7) 

Proo/. Let / G C^O^c), p be the solution of (1), (2), (3) and fix <r G S 1 . 
We define 

P(y,i) := / p(x,t)dL(x) 

and show that P is the unique solution of the two dimensional wave equation 
(4), (5), (6) with initial data F = (1l®l 3 )(f)(er, ■). Theorem 1 then follows 
from the fact that P(/)(<r, y rcc , t) = P(y rcc ,i) for y rcc G C rec and t > 0. 
From (2), (3) and the definitions of {TZ®Iz){f) and P it follows immediately 
that P satisfies the initial conditions (5), (6). It remains to prove that P 
satisfies (4). 

Since (ei(er), e2(<x), 63) is an orthonormal basis of M 3 for every x G M 3 
there exist unique elements s£R and y G K 2 such that x = sei(a) + y =: 
sei(cr) + yie2(<r) + ^2 e 3- Therefore, the Laplacian in M 3 decomposes into 
the sum 

92 A 

Using (1) results in 

f ( d 2 p d 2 p \ 

° = Jr V + y '^ ~ ds 2 ^ 1 ^ +y>*) ~ A y p(sei(cr) + y,t) 1 ds 

/ p(sei(er) + y,f)ds 



dt 2 



dp, , . . 
— (sei(o-) + y,t) 



- A y / p(sei(<r) + y,i) ds 
ii 

d 2 P 

= ^ 2 -(y,t)-A y P(y,t). 

The last equality holds since p(-, t) is compactly supported [12] for all t > 0. 
Hence P satisfies (4). □ 

The decomposition of V into a set of lower dimensional operators can be 
used to reduces the complexity of the derived reconstruction algorithms. 

The operator Q is closely related to the circular Radon transform with 
center set C rec [7, 18], which is defined by 

M : C °°(P» rec ) -> C°°(C rcc x (0,oo)) 

P ^ A* (P) := Uy rec , r)»^J F(y TCC + ru) dn(u)^j 
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dr (8) 



and maps a function F G (D TCC ) onto its integrals over circles with centers 
on the recording curve C rec - 

Using D'Alembert's formula the solution of (4), (5), (6) turns out to be 

[12] 

which means that Q(F) can be calculated if the circular Radon transform 
A4 (F) of F is known. The following Lemma states that also the reversal is 
true: 

Lemma 1. Let F G C$°(D rcc ). Then 

M(F)(y iec ,r) = -J q j=== dt (9) 

for all y rcc G C rcc and t > 0. 

A proof of a generalization of Lemma 1 is given in [23, Theorem 2.1.2], 
but for the sake of completeness we prove Lemma 1 since it is an elementary 
proof. 

Proof. Let F G C °°(i) rec ), fix y rcc G C rcc and define $(t) := -M(F)(y rec , t). 
Furthermore, we define for a function ^ G C°°((0, oo)) 

(J(*))(u) := f U -^=dt, 
Jo yu z — t z 

Jo Vu — t 
(£(*))(«) :=«(J(*))(«). 

Due to equation (8) it is sufficient to prove that (2/tt) J (( J7mui ($))') = 
<E>. As an auxiliary result we prove that for all G C°°((0, oo)) 

J( r ) : = (J (£(¥))) (r) = | ^ *(u)du . (10) 

For the proof of (10) let i > and consider the integral 

m = r * = r c f •m^wi ^ * . 



Jo \Jr 2 — t 2 Jo \Jo Vt 2 — u 2 J \Jr 2 — t 2 

From Fubini 's theorem we conclude that 
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We substitute t = s/v? + v 2 and dt = (v/t)dv in the inner integral. There- 
fore 

dv 



I(r) 



wo 



V r 2 — u 2 



\J r 2 — u 2 — v 2 



*5>(u)du 



Finally, we use the substitution v = \Jr 2 — u 2 sin (a) and dv = \Jr 2 — u 2 cos(a)da. 
We obtain 



I(r) 



r I /»tt/2 
\ JO 



7T 



da I $>(u)du = - I ^>{u)du 



Hence we have proved (10). 
Next we prove that 



[ID 



Using the relation 



^—\Jt 2 — r 2 
or 



and integration by parts we obtain 
r$(r) 



^—T 2 



dr 



$(r)Vt 



Vt 2 - r 2 



r=t 



+ / <$>'(r)\/t 2 - r 2 dr 



/ &(r)\/t 2 - r 2 dr . 
Jo 



'o 

By differentiating the above equation with respect to t we obtain 



Of 



+ t 



r=t Jo Vi 



dr . 



Hence we have proved (11). Applying J to (11) we find 

J '((,7mul(*))') = 3 = |*. 

For the last equality we applied the auxiliary result (10) to = <J>'. □ 

In order to prove that V is injective we use Lemma 1 and a uniqueness 
result for the circular Radon transform. 
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Theorem 2 (Uniqueness for arbitrary recording curve). Let V rec , D TCC and 
C rec be as above, and assume additionally, that D rec is convex and that C Tec C 
-D rec is a relatively open C 1 curve. If f G C^°(V rcc ) and V(f) = 0, then 



Proof. Due to (7) it is sufficient to show that 1Z and Q are injective. It is 
well known that the Radon transform 1Z is injective, see for example [17, 
Theorem 2]. Therefore, it remains to show that Q is injective. 

Let F G C£°(£> rec ) and Q(F) = 0. From (9) it follows that M(F) = 0. It 
is known that the circular Radon transform with center set D Tec is injective, 
see [2, Theorem 19] (compare also with [1, 7]). We conclude that F = 
which means that Q is injective. □ 

Equation (7) implies that thermoacoustic CT with line detectors deals 
with the inversion of the operator Q and the inversion of the two dimensional 
classical Radon transform 1Z. Fast and stable algorithms for inverting the 
two dimensional classical Radon transform 1Z are well investigated. The key 
task for deriving inversion algorithms for V is the inversion of the operator 
Q. Therefore, in the next Section we present and analyze an analytic formula 
for inverting Q. 

4 Linear recording curve: Inversion of the two di- 
mensional wave equation from data on a line 

Throughout the remaining parts of this article we consider the case where 
the recording curve C rec d M is a line. Depending on the angle between C rec 
and e2, the recording volume V rcc is either a closed cylinder, a closed cone 
or a half space, each considered as a subset of M 3 (see Figure 3). In all cases 
Q maps the initial data F onto the solution of of the two dimensional wave 
equation (4)- (6) restricted to C rec . 

From symmetry properties of the wave equation (bottom right image in 
Figure 3) it follows that in order to invert Q it is sufficient to solve the follow- 
ing problem: Recover H G Co°(]R x (0, oo)) from data G := -P|r x {o}x(o,oo)) 
where P is the solution of 



/ = 0. 




(12) 



11 




Figure 3: The array of line detectors measures the solution of the two di- 
mensional wave equation (4)- (6) restricted to C rec and is rotated in a plane 
(top left), tangential to a cylinder (top right) or tangential to a cone (bottom 
left). Bottom right: The restriction of the solution to C rec is uniquely deter- 
mined by the means over circles with centers on C Tec (see (8)) and remains 
unchanged after symmetrizing the initial data F around C rec . 

with symmetric initial conditions 

P(u,v,0) = H(u,\v\) , (u,v)£R 2 , (13) 
OP 

— (u,v,0)=0, (u,v)eR 2 . (14) 
u denotes the coordinate in direction C rec and v the signed distance to the 
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orthogonal direction. 

In the following we have to apply the Fourier-cosine transform to G. 
Since G ^ l'(R x (0, oo)), see Corollary 1 below, we have to sidestep to 
the space of symmetric tempered distributions. See Appendix A for the 
definitions of the spaces S sym (R x R) , C s ° y ° m (R x R) , S' sym (R x R) , the Fourier- 
cosine transforms C, C and the natural injections is, ic°° that map 
functions onto tempered distributions. We note that, since G sym G C°°(R x 
R), ic°° [Grgyni] G /Sg ym (Rxl) and the Fourier-cosine transform C'[ic°° [G S ym]] 
is well defined, see Appendix A for details. 

Theorem 3 (Fourier inversion formula). Let H G C^°(R x (0, oo)), P be 
the solution of (12)-(14), G := -P|r x {o}x(o,oo) an d ^ := G R 2 : oj > 

|A|}._ Then there exists G G L 1 (R x (0,oo)') n C°°(ft) with C'[i c °° [G sym ]] = 
lLipgym], and the Fourier inversion formula holds 



C[H](X, k)=G (A, a/^TA 2 ) • r^-^ (15) 
point-wise for all (A, k) G R x (0, oo). 

Proof. An elementary calculation shows that for all (A, k) G R x (0, oo) the 
functions 

cos(nv)e iXu cos (t^K 2 + A 2 



are solutions of (12) and (14). Since H G C^°(R x (0,oo)) C S' sym (R x R), 
H := C[H] is well defined. From the linearity of the wave equation and 
the inversion formula for the Fourier-cosine transform on 5 sym (R x R), see 
Proposition 3 in Appendix A, it follows that the unique solution of (12)-(14) 
is given by 

P(u, «, t) = 1 y H{\, as) cos (t \/k 2 + A 2 ) cos(kv) dn^j e iXu d\ . 

A 2 in the inner integral and taking v = leads to 

UJ 2 - A 2 , A) cos(Lut)^t= ) e iXu d\ (16) 
Vw 2 — A 2 / 



(17) 




H(X, v / w 2 -A 2 )cj/ v / ^ 2 -A 2 if u > |A| 
otherwise 
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G satisfies 



G\\ L r = J J™ \G{\,u)\ du dX = J QH \H (A, - A 2 ) l;^=p) ^ . 



Substituting w 2 = k 2 + A 2 in the inner integral we conclude that ||G||/,i = 
||i?|| L i and hence G G L X (R x (0,oo)). Equation (16) and G G L^M x 
(0,oo)) imply that C[G]~ = G and since the definitions of the Fourier- 
cosine transforms on S^^R x R) and L X (R x (0, 00)) are compatible, see 
Proposition 2 

[iLil^sym]] = ic°°[C[G]~ ym ] = i C °° [G sym ] . 

Applying the inversion formula for the Fourier-cosine transform on S'g ym (R x 
R), see Proposition 1, shows i L i [G sym ] = C'[ic°° [G sym ]]. 

Finally, solving (17) for H shows that (15) holds point-wise for all k = 
Vuj 2 - A 2 > 0. □ 

Corollary 1. Let H , G and G be as in Theorem 3. Moreover, assume that 
H 7^ is non-negative. Then G is not absolutely integrable. 

Proof. Assuming G G L 1 (R x (0, 00)) it follows from Proposition 3 in Ap- 
pendix A that G = C[G] G C°°(R x (0, 00)). Since H is non-vanishing it 
follows that C[i2](0, 0) > and since C[H] is continuous there exists e > 0, 
such that C[H](e,e) > 0. Inserting the sequence (e + 1/n, e + 2/n) in (17) 
shows that G is unbounded, contradicting G G C°°(R x (0,oo)). □ 

Corollary 1 implies that the integral 

- f [ G{u,t)cos{ujt)e~ iXu dtdu (18) 
n Jm. Jo 

is not absolutely convergent and therefore cannot be used directly to find 
an analytic representation of G. However, in the following we show that 
(18) exists in some generalized sense and can be used to find an analytic 
representation of G. 
In the following 

Fi[*](A) := (JA ' J^(u)e~ iXu du, A G R, 
/ 2 \ x / 2 r°° 

C 2 [$]H := I- J $(t)cos(ut)dt, u>0 
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denote the one dimensional Fourier transform on L 1 (M) and the Fourier- 
cosine transform on L 1 ((0, oo)), respectively. When applied to a function 
defined on M x (0, oo) Fi acts on the first and C2 on the second component. 
Finally, we recall the following special case of Jordan's theorem: If $ € 
L 1 ((0, 00)) is continuous differentiable on an interval including the point 
u > 0, then 

/o\V2 rK 

$(w)= lim - / cos(ut)C 2 M(t)dt . (19) 

K^oo\irJ J 

This enables us to prove our final result, namely an explicit expression 
for the calculation of G. 

Theorem 4 (Analytic formula for G). Let H, G, G and as in Theorem 
3. Then, for (A, to) G O 

G(A,w)= lim - f ( [ G{u,t)e- iXu du\ cos{ojt)dt . (20) 

K^oo TT Jq \J r J 

Proof From Theorem 3 it follows that Gei^Rx (0, 00)) and C[G] = G~ . 
Fubini's theorem implies that 

FioC 2 [G] = C[G] = G-. (21) 

Since the initial data H of the two dimensional wave equation (12)-(14) 
is compactly supported, G(-,t) £ Cg°(]R) for fixed t > and therefore 
Fi[G](-,t) G 5(M). Hence, applying Fx to (21) yields 

F X [G]- = C 2 [G]- (22) 

Next we apply (19) to G(A, •), for fixed A, 

G(X,u)= lim (-\ ' [ C 2 [G](\,t)cos(ut)dt . (23) 

K^oo \7Tj J 

In order to finish the proof we apply (22) to the integrand of (23) and 
obtain 

/ 2 \l/2 , K 

G(\,u) = lim - / Fi[G](A,t)cos(u;t)dt 
K->oo \irj J 

= lim - [ ( [ G{u,t)e~ iXu du\ cos(ut) dt . 

K^oo TT J \J R J 



□ 
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Theorem 4 shows that the calculation of the Fourier-cosine transform in 
the space S' m (R x R) of tempered distributions can be avoided. Moreover, 
in practical applications data G(u, t) can be collected on a finite space-time 
domain only. The fact that G(-,t) is compactly supported for t fixed and 
Theorem 4 justify that, in order to approximately calculate G, the domain 
of integration in (20) can be replaced by this finite domain. However, data 
truncation will introduce error in the reconstruction. 



INITIAL DATA H(u,v; 




PRESSURE G(u,t) 




20 40 60 80 100 
ct [mml 



RECONSTRUCTION 



RECONSTRUCTION 10 % noise 






Figure 4: Reconstruction with Fourier Formula. Top left. Initial 
data H. Top right. Simulated data with \u\ < 100 mm. Bottom left. 
Reconstruction from exact data. Bottom right. Reconstruction from noisy 
data. 



Remark 1 (Numerical realization). It is straight forward to derive a fast 
Fourier reconstruction algorithm based on Theorems 3, 4: First use the FFT 
algorithm to approximate G, see (20). Then use bilinear interpolation and 
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(15) to approximate C[F] on a cartesian grid. Finally, again by using the FFT 
algorithm, calculate an approximation of F. Due to the FFT algorithm, using 
data G to calculate F at N 2 grid points, the numer ical effort is 0(N 2 logN) 
only. 

Figure 4 shows results of a numerical example with initial data con- 
sisting of nine spheres supported in [—10, 10] x [0,20]. Data was collected 
with an array covering [—100, 100]. The solution of (12)-(14) was calculated 
analytically and uniformly distributed random noise with values in the in- 
terval [— e, e], where e equals to 5% of the largest data value, was added. 
The bottom images are obtained via the Fourier reconstruction formula, as 
presented in the Remark 1. 

The left image in Figure 5 shows a reconstruction where data was col- 
lected for u G [—20,20] only (limited aperture). In this case, boundaries of 
some spheres are not recovered correctly. We emphasize that this confirms 
theoretical predictions from microlocal analysis [16], namely that in order 
to stable recover a boundary with tangent t at Pj = (uj, Vi), data G(u, •) has 
to be collected such that (u{ — u,Vi) is orthogonal to t. As illustrated in 
Figure 5, this requirement is not satisfied for nearly horizontal boundaries. 




5 10 15 20 



v [mm] 

Figure 5: Finite aperture effect. Left. Reconstruction from data P(u, t) 
with \u\ < 20 mm. Nearly horizontal boundaries are not recovered correctly. 
Right. The boundary (Pi,ti) can be recovered stable, the boundary (P2,t2) 
cannot. 
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5 Conclusion and outlook 



In this paper we presented a mathematical framework of thermoacoustic CT 
with line detectors. The advantage of using line detectors is the improved 
imaging resolution compared to data acquisition based on point data. We 
proved that the three dimensional imaging problem reduces to the problem 
of reconstructing the initial data of the two dimensional wave equation from 
boundary measurements. 

In practical applications the output of a planar detector array offers 
only information from a finite aperture. Therefore, our future work will 
focus on the case of limited data. Furthermore, we will investigate scanning 
geometries where the array of line detectors forms different surfaces, e.g. 
cylinders which allow to collect a complete data set. 

A The Fourier— cosine transform on the space of 
symmetric tempered distributions 

In this appendix we define the Fourier-cosine transform on the space of sym- 
metric tempered distributions and prove the propositions stated in Section 
4. We will reduce the Fourier-cosine to the Fourier transform and therefore 
recall the basic definitions and properties of tempered distributions and the 
Fourier transform [21]. 

The Schwartz space S(M n x R) denotes the set of all infinitely differen- 
tiable functions ip : R n x R — > C that are rapidly decreasing together with 
all their derivatives, i.e. 

Vm £ No,/3 € Nq +1 : |M| m ,/3 := sup (1 + ||u'|r)|.D / Vu , )| < oo 

The variables in R n xR are denoted by u' = (u, v). We discriminate between 
u and v since the last component will play a special role. Equipped with the 
topology generated by the semi- norms ||v||m,/3j where m € No, (3 G Nq +1 , 
S(W l x R) is a Frechet space (metrizable and locally convex). 

The space S"(R n x R) of tempered distributions is defined as the topologi- 
cal dual of S, i.e. the set of all linear continuous functionals T : S(W l x R) — > 
C, equipped with the weak-* topology cr(S"(R n x R), 5(R n x R)). The eval- 
uation of a tempered distribution is denoted by (T,ip) := T(<p). The map- 
pings i s : S(R n x R) S"(R n x R), i L i : L 1 (R n x R) S"(R n x R), and 
i C oo : C°°(R n x R) -> S'{W l x R), where C°°(R n x R) denotes the set of 
all continuous functions tp with lrm|| u /i|_ KX3 y(u') = are well defined and 
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injective. In all cases a function tp is mapped to the distribution i[<p] defined 



where i stands either for ig, i L i or ic°°- 
Definition 1. 



• For tp : R n x (0, oo) -> C and ip : R n x R -> C define 



<£ sym : M" x 1 ^ C 

(u, u) i-> c^ sym (u, v) := ip(u, \v\) , 

(u, v) i ^ V- (u, u) := ~v) , 
f :R"xK^C 

(u, u) i-> if)- (u, u) := V(-u, v) , 
■01 : M" x 1 -» C 

(u, v) i — ► V- ( u ) u ) := VK — U > — u ) • 



• The Fourier transform on L 1 (IR n x R) is defined by 

F : L 1 (R n xi)-t C 00 (R" x R) 

/ / 1 \ (n+l)/2 , / , 



• The Fourier transform F : S(R n x R) -> ^(R" x R) restricted to the 
Schwartz space is well defined and a topological isomorphism (linear, 
bijective and bicontinuous) with inverse F _1 [<I>] = F[$]Z. 

• The Fourier transform on S"(R n x R) 



F' : S'(R n x R) -> S'(R n x R) 

T^F'[T] := (^<T,F[<d» 



is a topological isomorphism with ((F') 1 [T],(/?) = (T, F" 1 [y]). 
The definitions of the Fourier transform F and F' are compatible, i.e. for 



by 






F'[i L iM] =i C oo[FM] 



(24) 
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because, by using Fubini's theorem, 
(F / [i Ll M],^) = (i Ll M,F[ 



/ / 1 \ ("+l)/2 /■ \ 
p(x) (( — J J V(C)e^ <x '°(i n+1 C j d n+1 x 

= [ V>(0 ((f] in+1)/2 I ^(x)e- i ^OcP+ 1 x ] rf n+1 C 

= (i C oo[F^]]^} . 

In a similar way it can be shown that 

(FO-Hi Ll M]=ic°°[FM:]. (25) 

Definition 2. 

• A tempered distribution T G S"(R n x R) is called symmetric if 

V<^ G S(R n x R) : (T, = (T, if-) . 

• We define 

S Bym (R n x R) := {ip G C°°(R n x (0,oo),C) : ^ sym G S(R n x R)} , 
{ip G C°°(R n x (0,oo),C) : ^ sym G C°°(R n x R)} , 
{T G S'(R n x R,C) : T symmetric} . 



C s ° y ° m (R" x 

5'sym( IRn X 



The subspaces S' sym (R n xR), Sg (M n xR) have the topologies induced 
by S(R n x R) and S'(R n x R). 

• The Fourier-cosine transform C on L 1 (R n x (0,oo)) is defined by 

/ i \ (n+l)/2 p / poo \ 

C[ip](X, k) := 2 i—J J Ij <p(u, v) cos(Kv)dv) e~ l{X ' u) d n u , 

• The Fourier-cosine transform C := F' |^ m (R»xR) on S' syra (R n x R) 
denotes the restriction of F' : S'(R n x R) -> 5"(R n x R) to the space 
S sym (R"xR). 

Proposition 1 (C is an isomorphism). The mapping C : S' syni (R n x R) — > 
Sg ym (l" x R) is a topological isomorphism and ((C)' 1 [T], ip) = (C[T],(p~). 
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Proof. The mapping C : S' sym (R n xl)-t S'(R n x R) is linear, continuous 
and injective. Let T U T 2 G S^ ym (R n x R) and v? G S(R n x R). Then 

(F'tTi],^} = (ri,F[y>_]> = (Tx.FM-) = (T 1? FM) = (F'pU^) 
((F , )- 1 [T 2 ],^_> = (r 2 ,F- 1 [^_]} = (r 2 ,F- 1 M_) = <T 2 ,F>]} = ((F / )- 1 [T 2 ],^} 

Hence, F'[Ti] and (F')" 1 ^] are symmetric and C : S^ ym (R n xR) ^ ym (R n x 
R) is well defined, surjective and therefore an isomorphism. Finally, 

{{C')- l [TiU) = m-'lTi}^) = (Ti, F" 1 [tp]) = (Tx,F[ipZ\) 

= (Ti,F[(pZ]-) = (Ti,F[(p~]) = (F'in},^) = (C[ri]^-> 

finishes the proof. □ 

Lemma 2. If ip G L 1 (R" x (0,oo)), </> G C° y ° m (R n x K), v G 5 sym (R n x R) 
then i L i^sym],ic oo [^ sy m],is[?/sym] G 5^ ym (R n x R) and F[tp syni ] = C[ip] sym . 

Proof. Let 7 G tp, rf\ and i denote the corresponding embedding. Then 

(i[7sym],p) = / 7 sym (u')p(u / )(i n u / 

7(u, |v|)p(u, f )dv I <i n u 
/ 

7(u, \v\)p(u, -t>)<2u ) cTu = (i[7 sym ],p_) . 



which proves the first statement. Moreover, for k > 

(n+l)/2 



C[<p](\,K) = 2 f — J y M ^(u,u)cos(Ku)duJ e- i(A ' u) d n u 

/ ( f <p syra (u,v)e- iKV dv) e-^cTu 

JE" UR / 



1 N (n+l)/2 

= F[^ sym ](A, k) . 

ip sym implies that F[(^ sym ] is a symmetric function with respect to the last 
component, i.e. F[y? sym ] G C°°(R n x R) symmetric. So, by extending C[(p] 
to C[y] sym the second statement holds. □ 

Proposition 2 (Compatibility). If ip G L 1 (R n x (0, 00)) then C[tp] sym G 
C s ° y ° m (R"xR). Moreover, C'[i L ifoym]] = i C oo[C[^] sym ] and (C)- 1 [i L i [ip sym ]) = 
ic°°[C[<p]- m ]. 
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Proof. Since <^ sym G L x (R n x M) and F'[i L i [<^ sym ]] = i c °o [F[^ sym ]] (see (24)) 
the restriction of F' to the space S' syra yields 

CpH^gym]] = ic°°[ F bsym]] • (26) 

Furthermore, by using (25) 

(C'^fiLiMsym] = ic°°[F[</? sym ]l] = ic°°[F[v3 S ym]~] • (27) 

Applying Lemma 2 to (26) and (27) accomplishes the proof. □ 

Proposition 3 (C is an isomorphism). The mapping C : /S sym (R n x R) — > 
<S' sym (M n x R) is a topological isomorphism with 

C- 1 [tp](u,v) = 2(^\ V j (J tp(\, k) c0s(Kv)di?J e i{x ' u) d n \ . 

Proof. Both statements follow from C[</?] sym = F[<^ sym ]| R n x (o >OQ \ and the 
corresponding result for F. □ 

We finally note that analogous constructions are possible for the space 
5(R n x R' m ) of all Schwartz functions that are rotationally invariant in the 
second component. This leads to an extension of the Fourier-Hankel trans- 
form to a space of distributions. In the special situation m = n+1 this space 
was denoted by 5^(R n x R n+1 ) in [3, 22] and used to invert the spherical 
mean operator. 

B Inversion of the wave equation in arbitrary di- 
mensions. 

Let n G N, n > 1, and H G C§°(R n x (0, oo)). The variables in R n x R are 
once more denoted by u' = (u, v). We consider the problem of reconstructing 
the initial data H of the (n + l)-dimensional wave equation 

d 2 \ 

— -A u , P = 0, (u, v, t) G R n x R x (0, oo) (28) 
at J 

with initial conditions 



P(u',0) =H(u,\v\), u'ePxl, (29) 
dP 

— (u ', 0) = , u' G R n x R (30) 
from data G(u,t) := P(u,v = 0,t). Here A u / denotes the Laplacian in 
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Theorem 5. Let H G Cg°(R n x (0,oo)), P the solution of (28)-(30) and 
G := PI 

R n x{o}x(o,oo) denote the restriction of P to W l x {0} x (0, oo). 
Then there exists G £ L 1 (R n x (0, oo)), such that C'[ic°° [G sym ]] = 
i L i[G sym \. Moreover, 



f 1 \ (n+l)/2 ,K / p \ 

G(X,lo)= lim 2 — / / G(u, t)e- l ^d n u ) cos(wt) dt 

K->oo \2irJ J Q \J R n ) 

exists point-wise for all uj > A and 

C[H] (A, K ) = G (A, v^ 2 + I|A|| 2 ) • . 9 K „ xll9 , « > . 
v 7 y k z + ||A||^ 

The proof of Theorem 5 is analogous to the corresponding result in n + 
1 = 2 dimensions, as presented in Section 4. 
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